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Abstract 

The Quantitative Rescattering Theory (QRS) for high-order harmonic generation (HHG) by 
intense laser pulses is presented. According to the QRS, HHG spectra can be expressed as a prod- 
uct of a returning electron wave packet and the photo-recombination differential cross section of 
the laser-free continuum electron back to the initial bound state. We show that the shape of the 
returning electron wave packet is determined mostly by the laser only. The returning electron 
wave packets can be obtained from the strong-field approximation or from the solution of the time- 
dependent Schrodinger equation (TDSE) for a reference atom. The validity of the QRS is carefully 
examined by checking against accurate results for both harmonic magnitude and phase from the 
solution of the TDSE for atomic targets within the single active electron approximation. Com- 
bining with accurate transition dipoles obtained from state-of-the-art molecular photoionization 
calculations, we further show that available experimental measurements for HHG from partially 
aligned molecules can be explained by the QRS. Our results show that quantitative description of 
the HHG from aligned molecules has become possible. Since infrared lasers of pulse durations of a 
few femtoseconds are easily available in the laboratory, they may be used for dynamic imaging of 
a transient molecule with femtosecond temporal resolutions. 

PACS numbers: 42.65.Ky, 33.80.Rv 
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I. INTRODUCTION 



High-order harmonic generation (HHG) has been studied extensively since 1990's, both 
experimentally and theoretically. Initial interest in HHG was related to the generation of 
coherent soft X-ray beams Q, 2], which are currently being used for many applications in 



ultrafast science experiments [3|, |4|. In the past decade, HHG has also been used for the 



production of single attosecond pulses 



and attosecond pulse trains 



thus opening 



up new opportunities for attosecond time-resolved spectroscopy. 

While HHG has been well studied for atoms, much less has been done for molecules. Initial 
interest in HHG from molecules was due to the fact that it offers promising prospects for in- 
creasing the low conversion efficiency for harmonic generation. The presence of the additional 
degrees of freedom such as the alignment and vibration opens possibilities of controlling the 
phase of the nonlinear polarization of the medium and of meeting the phase-matching condi- 
tion. Thanks to the recent advancements in molecular alignment and orientation techniques 
91], investigation of the dependence of HHG on molecular alignment reveals the distinctive 



features of molecular HHG [1(1 as well as their structure 



ll| . In particular, the existence of 



distinctive minima in the HHG spectra from H^ have been theoretically predicted by Lein et 



al. 
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ISj. The experimental measurements within the pump-probe scheme for CO2 indeed 



showed the signature of these minima for partially aligned CO 
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151]. More recently. 



newer experiments by JILA, Saclay, and Riken groups [161, Il7|, Il8|, [19] began to focus on 
the phase of the harmonics, using mixed gases and interferometry techniques. These studies 
have revealed that near the harmonic yield minimum the phase undergoes a big change. In 
view of these experimental developments, a quantitative theory for HHG from molecular 
targets is timely. 

HHG can be understood using the three-step model 
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2l| . First the electron is released 



into continuum by tunnel ionization; second, it is accelerated by the oscillating electric field 
of the laser and later driven back to the target ion; and third, the electron recombines with 
the ion to emit a high energy photon. A semiclassical formulation of this three-step model 
based on the strong field approximation (SFA) is given by Lewenstein et al 22|. In SEA, 
the liberated continuum electron experiences the full effect from the laser field, but not from 
the ion. In spite of this limitation, the model has been widely used for understanding the 
HHG by atoms and molecules. Since the continuum electron needs to come back to revisit 
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the parent ion in order to emit radiation, the neglect of electron-ion interaction is clearly 
questionable. Thus, over the years efforts have been made to improve upon the SFA model, 

2^ 



by including Coulomb distortion 
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251 ] ■ These improvements, when applied to simple 



systems, however, still do not lead to satisfactory agreement with accurate calculations, 
based on numerical solution of the time-dependent Schrodinger equation (TDSE). 

With moderate effort, direct numerical solution of the TDSE for atomic targets can be 
carried out, at least within the single active electron (SAE) approximation. For molecules. 



accurate numerical solution of the TDSE is much more computationa 



not been carried out except for the simplest molecules such as H2' 
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demanding and has 
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27l |. Thus most 



of the exis ting calculations for HHG from molecules were performed using the SFA model 



28|, 



29 



30 



31 



32 



3, 



34 



|. Additionally, in order to compare single molecule calculations 
with experimental measurements, macroscopic propagation of the emitted radiation field in 
a gas jet or chamber needs to be carried out, with input from calculations involving hundreds 
of laser intensities to account for intensity variation near the laser focus. For this purpose 
the TDSE method is clearly too time-consuming even for atomic targets. 

In view of the inaccuracies of the SFA and the practical inefficiency of the TDSE method. 



we have proposed the quantitative rescattering t 



wave SFA or SW-SFA in our earlier papers 



leory or the QRS (also called the scattering- 



36 



37( 1) as a simple and practical method for 



obtaining accurate HHG spectra. The main idea of the QRS is to employ exact transition 
dipoles with scattering wave instead of the commonly used plane-waves for the recombina- 
tion step in a HHG process. Within the QRS, induced dipole moment by the laser can be 
represented as a product of the returning electron wave packet and the complex recombina- 
tion transition dipole between free electrons with the atomic or molecular ion 
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. The 



shape of the wave packet has been shown to be largely independent of the target. Therefore 
the returning wave packet can be calculated from the standard SFA model or by solving 
the TDSE from a reference atom with similar ionization potential as of the molecule under 
consideration. The QRS approach has been shown to be much more accurate than the stan- 

n n 

dard SFA model for rare-gas atoms 361] and a prototypical molecular system H^ (371]. We 



note that the QRS has also been success: 
from various atomic systems 
of atoms. 
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40 



41 



uUy used to explain high-energy electron spectra 
4^ 



431 ]. and non- sequential double ionization [44] 



In order to use the QRS, we need to know the complex transition dipoles from aligned 
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molecules. In other words, we need both the differential photoionization (or photo- 
recombination) cross section and the phase of the transition dipole for a fixed-in-space 
molecule. Calculations for these quantities can be carried out using the available molecular 
photoionization methods, which have been developed over the last few decades. In this 
paper we employ the iterative Schwinger variational method developed by Lucchese and 
collaborators 
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46|. 



The goal of this paper is to give a detailed description of all the ingredients of the QRS 
and demonstrate its validity. As examples of its application we consider HHG from aligned 
O2 and CO2 and show that the QRS is capable of explaining recent experiments. The rest of 
this paper is arranged as follows. In Sec. II we summarize all the theoretical tools needed for 
calculating HHG spectra from partially aligned molecules. We also explain how photoion- 
ization (photo-recombination) differential cross sections and phases of the transitions dipole 
are calculated for both atomic and molecular targets. The QRS model is formulated and 
its validity is carefully examined in Sec. HI for atomic targets where results from the TDSE 
serve as "experimental" data. In Sec. IV we illustrate the application of the QRS method 
by considering the two examples of O2 and CO2 and compare the QRS results with available 
experiments. Sec. V discusses the predictions of the QRS for the minima in HHG spectra in 
relation to the simple two-center interference model. We also extend the QRS calculation 
for O2 to higher photon energies by using 1600-nm lasers, to explore the interference effect 
in HHG spectra from O2. Finally we finish our paper with a summary and outlook. 

We note that the macroscopic propagation effect has not been treated in this paper . 
In the literature, theoretical treatments of such effects have been limited to atomic targets, 
and mostly starting with the SEA to calculate single atom response to the lasers. In this 
respect, the QRS offers an attractive alternative as the starting point since its calculation 
of single atom response is nearly as fast as the SFA, but with an accuracy much closer to 
the TDSE. Atomic units are used throughout this paper unless otherwise stated. 



II. THEORETICAL BACKGROUND 



The theory part is separated into six subsections. We will describe all the ingredients used 
in the QRS. Eirst we consider the TDSE and the SEA, as the methods to calculate HHG 
spectra and to extract returning electron wave packets. We will then present theoretical 
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methods for calculation of photoionization (photo-recombination) cross sections for both 
atomic and linear molecular targets. We will also briefly describe intense laser ionization, as 
ionization rates are used in the QRS to account for depletion effect and overall normalization 
of the wave packets. Lastly, we explain how the partial alignment of molecules by an aligning 
laser is treated. 



A. Method of solving the time-dependent Schrodinger equation 



The method for solving the time-dependent Schrodinger equation ( 
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DS 
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) for an atom 



491]. Here we 



in an intense laser pulse has been described in our previous works 
present only the essential steps of the calculations and the modifications needed to treat 
the HHG problem. At present accurate numerical solution of the TDSE for molecules is 
still a formidable cornputational task and has been carried out mostly for simple molecular 



systems such as 12, l26|, l27[ . 



We treat the target atom in the single active electron model. The Hamiltonian for such 
an atom in the presence of a linearly polarized laser pulse can be written as 

H = Ho + H,{t) = -iv^ + V{r) + H,{t). (1) 

The atomic model potential V{r) is parameterized in the form 

Vir) = , (2) 

which can also be written as the sum of a short-range potential and an attractive Coulomb 
potential 

V{r) = ~- + Vs{r), (3) 
r 

The parameters in Eq. ([2]) are obtained by fitting the calculated binding energies from this 
potential to the experimental binding energies of the ground state and the first few excited 
states of the target atom. The parameters for the targets used in this paper can be found 
in 5^ . One can also use a scaled hydrogen as a reference atom, as shown in Ref . 0, [s^] ■ 
Here the effective nuclear charge is rescaled such that the ionization potential of the rescaled 
H(ls) matches the Ip of the atom or molecule under consideration. 
The electron-field interaction Hi{t), in length gauge, is given by 

H,{t) = r-E{t). (4) 



For a linearly polarized laser pulse (along the z axis) with carrier frequency lo and carrier- 
envelope-phase (CEP), the electric field is taken to have the form 

E{t) = zEq cos^ (—) cos{ujt + if) (5) 



r 

for the time interval (— r/2, r/2) and zero elsewhere. The pulse duration, defined as the full 
width at half maximum (FWHM) of the intensity, is given by F = r/2. 75. 

The time evolution of the electronic wavefunction \l/(r,t), which satisfies the TDSE, 

z^^^{r,t) = H^{r,t) (6) 

is solved by expanding in terms of eigenfunctions, Rni{r)Yim{r) , of Hq, within the box of 

r e [0,rmax] 

^{r,t) = J2cni{t)Rni{r)YUr) (7) 

nl 

where radial functions Rni{r) are expanded in a discrete variable representation (DVR) 
basis set [51] associated with Legendre polynomials, while Cni are calculated using the split- 
operator method [0 



n'l' 



xe-'^'^^'^'y^iy.cMt)- (8) 

For rare gas atoms, which have p-wave ground state wavefunctions, only m = is taken into 
account since for linearly polarized laser pulses, contribution to the ionization probability 
from m = ±1 is much smaller in comparison to the m = component. 

Once the time-dependent wavefunction is determined, one can calculate induced dipole 
in either the length or acceleration forms 

DL{t) = {^{r,t)\z\^{r,t)) (9) 

D^{t) = {^{r,t)\^^\^{r,t)) (10) 

For low intensities, when the ionization is insignificant, the two forms agree very well. For 
higher intensities the above length form should be corrected by a boundary term, to account 
for the non- negligible amount of electron escape to infinity, see Burnett et al [53] . Therefore, 
we found it more convenient to use the acceleration form. To avoid artificial refiection due to 
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a finite box-size, we use an absorber of tlie form cos^ ^ '^[k {r — r cut) I'^ij' max — ^cut)] for r > Vcut 
52], to filter out tlie wave packet reacliing tlie boundary. We typically use r^ax = 200 to 
400, Tcut ~ fmax — 100, witli up to about 800 DVR points and 80 partial-waves. We liave 
cliecked tliat tlie results are quite insensitive with respect to the absorber parameter r^ut- 



B. The strong field approximation (SFA) 



The strong- field approximation (SFA) has been widely used for theoretical simulations of 



HHG from atoms 22] and molecules 28 



29 



30 



31 



32 



33 



34 



35| . This model is known 



to give qualitatively good results, especially for harmonics near the cutoff. However, in the 
lower plateau region the SFA model is not accurate 36|]. Nevertheless, since the propagation 
of electrons after tunnel ionization is dominated by the laser field, one can use the SFA to 
extract quite accurate returning electron wave packet which can be used in the QRS theory. 
Here we briefly describe the SFA, extended for molecular targets 281 ]. 

Without loss of generahty, we assume that the molecules are aligned along the x-axis, in 
a laser field E{t), linearly polarized on the x-y plane with an angle 9 with respect to the 
molecular axis. The parallel component of the induced dipole moment can be written in the 
form 

/•oo / 7J- \ 



x[cos6d:rit - t) + smedy{t - T)]E{t - t) 
X exp[—iSst(t,T)]a*(t)a(t — r) + c.c. 



(11) 



where d{t) = d[pst{t, r)+A(t)], d{t—T) = d[pst{t, r)+A(t— r)] are the transition dipole mo- 
ments between the ground state and the continuum state, and Pst{t,T) = — A{t')dt' / t 
is the canonical momentum at the stationary points, with A the vector potential. The 
perpendicular component Dj_(t) is given by a similar formula with [cos 9dl(t) + sm9d*{t)] 
replaced by [sin 6*^* (t) — cos6d*(t)] in Eq. (ITT!) . The action at the stationary points for the 
electron propagating in the laser field is 

£ nPst{t,T)+A{t')]' ■ 



Sst{t,r) 



+ h dt', 



(12) 



where Ip is the ionization potential of the molecule. In Eq. (fTTj) . a{t) is introduced to account 
for the ground state depletion. 
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The HHG power spectra are obtained from Fourier components of the induced dipole 
moment D{t) as given by 



P(uj) oc \a(uj) 



(13) 



In our calculations we use ground state electronic wavefunctions obtained from the gen- 
eral quantum chemistry codes such as Gamess ^] or Gaussian [55|. Within the single active 
electron (SAE) approximation, we take the highest occupied molecular orbital (HOMO) for 
the "ground state". In the SFA the transition dipole d{k) is given as (fc|r|\E'o) with the con- 
tinuum state approximated by a plane-wave \k). In order to account for the depletion of the 
ground state, we approximate the ground state amplitude by a(t) = exp[— /l^ W{t') /2dt'], 
with the ionization rate W{t') obtained from the MO-ADK theory or the MO-SFA [see, 
Sec. 11(E)]. 



C. Calculation of transition dipoles for atomic targets 

Photo-recombination process is responsible for the last step in the three-step model for 
HHG. The QRS model goes beyond the standard plane- wave approximation (PWA) and 
employs exact scattering wave to calculate the transition dipoles. Since photo-recombination 
is the time-reversed process of photoionization, in this and the next subsection we will 
analyze the basic formulations and methods for calculating both processes in atoms and 
linear molecules. 

The photoionization cross section for transition from an initial bound state to the final 
continuum state \E'^ due to a linearly polarized light field is proportional to the modulus 
square of the transition dipole (in the length form) 

dk,r^iuj) = {^^^\r■n\<^'^). (14) 

Here n is the direction of the light polarization and k is the momentum of the ejected 
photoelectron. The photoionization differential cross section (DCS) can be expressed in the 



general form as 



5f 



{^,\vn\^-^) , (15) 



dfl^d^n c 

where + Ip = u with Ip being the ionization potential, u the photon energy, and c the 
speed of light. 



To be consistent with the treatment of the TDSE for atoms in laser fields in Sec. 11(A), 
we will use the model potential approach for atomic targets. The continuum wavefunction 
\l/;^(r) then satisfies the Schrodinger equation 



2 ^ ^ 2 



^'^{r) = 0, (16) 



where the spherically symmetric model potential V{r) for each target is the same as in 
Eq. m- 

The incoming scattering wave can be expanded in terms of partial waves as 

oo I 

= -TtT. E ^'exp[-^{al + 5l)]REl{r)Y^minr)YrJnk). (17) 

V 1=0 m=-l 

Here, Si is the l-th partial wave phase shift due to the short range potential Vs{r) in Eq. (I2l), 
and ai is the Coulomb phase shift 

ai = argr(/ + l + n), (18) 
7 = ~Z/k, (19) 

with the asymptotic nuclear charge Z = 1. Rei is the energy normalized radial wavefunction 
such that 

roc 

/ REiir)REnirydr = 5{E-E'), (20) 
Jo 

and has the asymptotic form 

REi{r) ^ i^sin(A;r - l7i/2 --f\og2kr + ai + 6i). (21) 
The initial bound state can be written as 

^,ir)=R^iXr)Yi,^A^r), (22) 
\Rr,ifr^dr = l. (23) 
In our calculations we solve both the bound state and scattering states numerically to obtain 



Rnii and Rei- In the PWA, the continuum electron is given by the plane waves 

^fc(0 = ^^^exp(«fc-r). (24) 

where the interaction between the continuum electron and the target ion in Eq. (fT6|) is 
completely neglected. 
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It is appropriate to make additional comments on the use of PWA for describing the con- 
tinuum electron since it is used to calculate the dipole matrix elements in the SFA. Within 
this model, the target structure enters only through the initial ground state. Thanks to 
this approximation, the transition dipole moment is then given by the Fourier transform of 
the ground state wavefunction weighted by the dipole operator. Thus by performing inverse 
Fourier transform, the ground state molecular wavefunction can be reconstructed from the 
transition dipole moments. This forms the theoretical foundation of the tomographic pro- 



cedure used by Itatani et al However, it is well-known that plane wave is not accurate 
for describing continuum electrons at low energies (~ 20 eV up to ~ 0.5 keV), that are 
typical for most HHG experiments, and that all major features of molecular photoionization 
in this energy range are attributable to the property of continuum wavefunction instead of 
the ground state wavefunction. In fact, the use of PWA in the SFA model is the major 
deficiency of the SFA which has led to inaccuracies in the HHG spectra, as shown earlier in 



36 



371] and will also be discussed in Sec. HI. 
Note that so far we have considered one-photon photoionization process only. As mention 
earlier, the more relevant quantity to the HHG process is its time-reversed one-photon photo- 
recombination process. The photo-recombination DCS can be written as 

__ = _|(*.|,.„|<p.)| . (25) 

In comparison with photoionization DCS in Eq. fllSI) . apart from a different overall factor, 
the continuum state here is taken as the outgoing scattering wave instead of an incoming 
wave \E'^. The partial wave expansion for \I'^ is 

1 oo / 

Hir) = ^T.T. ^'exp[^ial + 5l)]RElir)Y^minr.)Y:^ink). (26) 

V 1=0 m=-l 

We note here the only difference is in the sign of the phase (cr^ + 6i) as compared to Eq. (1171) . 
In fact, the photoionization and photo- recombination DCS's are related by 

= (27) 
which follows the principle of detailed balancing for the direct and time-reversed processes 



57|. 



For clarity in the following we discuss the photo-recombination in argon. Once the scat- 
tering wave is available, the transition dipole can be calculated in the partial- wave expansion 

11 



as 



/ i{o-;+5i) 



{Rnl^\r\REl) 



Im 

x{Yi^^^\cos9\Yi^)Y:jn,,). 



(28) 



Here the polarization direction n is assumed to be parallel to z-axis. Using the relation 

1^ 



cos 6* 



the angular integration can be written as 




Y 



101 



(29) 



{Yi^m^l cos 9\Yi^) = ^J YC^^{9,<p)Y,o{9,mm{0A)siri9d9d^, (30) 



which can also be expressed in terms of the Wigner 3j-symbol by using 

Yi,^,i( 



)Y,, 



sm 



'(2/i + 1) (2/2 + 1) (2/3 + 1) 



X 



h h h 



An 
I 



h h h 
mi 1712 fns 



(31) 





From this equation, it is clear that only m = rrii and I = k — 1 and I = k + 1 contribute. 
This gives [58 1 



A'^m. = {Yk^.\ cos 9\Yi^^) 



{2li+l){2h-l) 



for I = h — 1 



for / = /i + 1 



(32) 



(2Zi+3)(2ii+l) 

As discussed in Sec. 11(A), we only need to consider the case of recombination of electron 
back to Ar(3po), as electrons in m = ±1 states are not removed by tunnel ionization in the 
first step. From the above equation we see that only / = /j — 1 = (s-wave) and I = U + l = 2 
(d-wave) contribute, with = l/Vs and Af^ = 2/vT5, respectively. Furthermore, most 
contribution to the HHG process comes from electrons moving along laser's polarization 
direction. Since 

WO) = /J, 

>^;o(o,o) = /J, 



for A; II n the transition dipole can be written as 

1 



i((To+5o) 



{R3i\r\REo)/2 



-e*('^^+'^^)(i?3i|r|i?s2) 



(33) 
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Note that the transition dipole is intrinsically a complex number. The dominant com- 
ponent is the d-wave. Thus when the real d-wave radial dipole matrix element vanishes, 

n 

the cross section will show a minimum known as the Cooper minimum [59]. Note that the 
cross section does not go precisely to zero because of the contribution from s-wave (the first 
term). We will analyze this example in more detail in Sec. III(B). 

We comment that the calculation of transition dipole moment for atomic targets presented 
above is based on the SAE. The validity of such a model for descri bing photoionization has 



been studied in the late 1960's, see the review by Fano and Cooper j60|. Such a model gives 
an adequate description of the global energy dependence of photoionization cross sections. 
To interpret precise photoionization cross sections such as those carried out with synchrotron 
radiation, advanced theoretical methods such as many-body perturbation theory or R-matrix 
methods are needed. For HHG, the returning electrons have a broad energy distribution as 
opposed to the nearly monochromatic light from synchrotron radiation light sources, thus 
the simple SAE can be used to study the global HHG spectra. 

D. Calculation of transition dipoles for linear molecules 

Let us now consider photoionization of a linear molecule. For molecules, the complication 
comes from the fact that the spherical symmetry is lost and additional degrees of freedom 
are introduced. The photoionization DCS in the body-fixed frame can be expressed in the 



general form 45 1 



(Pa Air'^uk , , , , ,2 /„ 

|4,r^(^)| • (34) 



For the emitted HHG component with polarization parallel to that of the driving field, 
the only case to be considered is || n. To treat the dependence of the cross section on 
the target alignment, it is convenient to expand the transition dipole in terms of spherical 
harmonics 

/47r\^/2 

d^A^) = f y j E dlmMYl*mi^>^)yi,i^r.). (35) 

Here the partial-wave transition dipole is given by 

dim>.{^) = {^i\r^\^kiJ, (36) 
with ru = z for linear polarization. 
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In our calculations, we use an initial bound state obtained from the MOLPRO code 



6l| within the valence complete-active-space self-consistent field (VCASSCF) method. The 
final state is then described in a single-channel approximation where the target part of the 
wavefunction is given by a valence complete active space configuration interaction (VCASCI) 
wavefunction, obtained using the same bound orbitals as are used in the wavefunction of 
the initial state. The Schrodinger equation for the continuum electron is 



+ V(r) 

2 r ^ ' 2 



r 



0, (37) 



where V{r) is the short-range part of the electron molecule interaction, which will be dis- 
cussed below. Note that the potential is not spherically symmetric for molecular systems. 
The Schrodinger equation fl57|) is then solved by using the iterative Schwinger variational 
method [45|]. The continuum wavefunction is expanded in terms of partial waves as 

<Pkir) = (lY'i: E UMmir)Y;jnk), (38) 

l=Om=-l 

where an infinite sum over / has been truncated at / = Ip. In our calculations, we typically 
choose Ip = 11. Note that our continuum wavefunction is constructed to be orthogonal to 
the strongly occupied orbitals. This avoids the spurious singularities which can occur when 



scattering from correlated targets is considered 62|]. We have used a single-center expansion 
approach to evaluate all required matrix elements. That means that all functions, including 
the scattering wavefunction, occupied orbitals, and potential are expanded about a common 
origin, taken to be the center of mass of the molecule, as a sum of spherical harmonics times 
radial functions 

E E fim{r)YiUOA)- (39) 

1=0 m=-l 

With this expansion, the angular integration can be done analytically and all three- 
dimensional integrals reduce to a sum of radial integrals, which are computed on a radial 
grid. Typically, we use Imax = 60 to 85. 

Next we describe how the interaction potential, V is constructed. The electronic part of 
the Hamiltonian can be written as 

N N 1 

^ = EMO + E- (40) 

with 



i=l i<j 



MO = -f-E^, (41) 

^ rt ' in. 
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where Za are the nuclear charges and is the number of electrons. In the single channel 
approximation used here, the ionized state wavefunction \E'fc is of the form 



^fc = A($0fc) (42) 

where $ is the correlated — 1 electron ion core wavefunction, 0^ is the one-electron con- 
tinuum wavefunction, and the operator A performs the appropriate antisymmetrization and 
spin and spatial symmetry adaptation of the product of the ion and continuum wavefunc- 
tions. The single-particle equation for the continuum electron is obtained from 

{5^fc|i/-E|^fc) = 0, (43) 

where b'^k is written as in Eq. fl42p . with 0^ replaced by (50fc- By requiring this equation 



to be satisfied for all possible 5\E'fc (or 50fc), one obtains a non-local optical 



63, 



potential which 



6J] as written in 



can be written in the form of a Phillips-Kleinman pseudopotential, V ^ 
Eq. dSZD 

Note that ionization from a molecular orbital other than HOMO, say HOMO-1, can be 
done in the same manner, except that the target state $ employed in Eq. (l42l) needs to 
be replaced by the wavefunction for the corresponding excited ion state the corresponds 
to ionization from the HOMO-1 orbital. Furthermore, the above single-channel formalism 
can be extended to coupled-multichannel calculations to account for additional electron 



correlation effects 



461]. In this paper we limit ourselves to single-channel calculations. 



This single-center expansion approach has also been implemented for non-linear targets 



in the 
tia 



rozen-core Hartree-Fock approximation including the full non-local exchange poten- 



65l |. A somewhat simpler approach, the finite-element R-matrix (FERM3D) by Tonzani 



66| can also be employed to calculate transition dipoles and photoionization cross sections. 



The FERM3D code is especially well adapted for complex molecules. In FERM3D the elec- 
trostatic potential is typically obtained from general ah initio quantum chemistry software 
such as GAUSSIAN ^] or GAMESS 5J] and the exchange potential is approximated using 
a local density functional. A polarization potential is also added to describe the long range 
attraction between the continuum electron and the target ion. This code can also calculate 
ionization from any occupied molecular orbitals. 
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E. Description of strong-field ionization from the MO-ADK and the MO-SFA 



In the tunneling regime, the most successful general theories for ionization from molecules 



are the MO-ADK 



67l | and MO-SFA [68|, l69|. For our purpose of studying HHG process, 



it is important to make clear distinction between the total ionization (integrated over all 
emission directions) and the (differential) ionization along the laser polarization direction. 
The former is used in the description of the ground state depletion [see Sec. 11(B)], while the 
latter is directly related to the magnitude of the returning electron wave packet, which will 
be described in Sec. Ill (A). It is well-known that the SFA (or MO-SFA) can give qualitatively 
good ATI spectra, but not the overall magnitude |48|, whereas the MO-ADK can give total 
ionization rates. Strictly speaking, the MO-ADK only gives total ionization, whereas the 
MO-SFA can also give differential rates. Therefore in our calculations we use both theories. 
For N2 and O2, the total ionization rates from the MO-ADK and the (renormalized) MO-SFA 
agree very well 7^. However for CO2 the alignment dependent rates from the MO-ADK, 
the MO-SFA, and the recent experiment TjJ all disagree with each other, with the MO-ADK 
predicting a peak near 30°, conapared to 40° from the MO-SFA theory, and a very sharp 
peak near 45° from experiment 71 1. 



The MO-ADK theory is described in details 67|]. Here we will only mention briefly main 
equations in the MO-SFA theory. In the MO-SFA model the ionization amplitude for a 
transition from a bound state $o(^) to continuum is given by 22] 

/oo 
dt{p + A{t)\r ■ E{t)\^o) exp[-tS{p, t)], (44) 
-00 

where 

S(p.*)=f + (45) 

with p the momentum of the emitted electron, Ip the binding energy of the initial state, and 

A{t) the vector potential. In the SFA the effect of the core potential is totally neglected in 

the continuum state, which is approximated by a Volkov state 

1 



{r\p + A{t)) 



■ exp {i [p + A{t)] ■ r} . 



(46) 



(27r)3/2 

For the bound states we use the wavefunctions generated from the ab initio quantum chem- 
istry Gaussian [551] or Gamess codes (54!]. In this regard we note that the use of these 
wavefunctions in the MO-ADK might not be sufficiently accurate, as the MO-ADK needs 
accurate wavefunctions in the asymptotic region. 
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F. Alignment distributions of molecules in laser fields 



When a molecule is placed in a short laser field (the pump), the laser will excite a 
rotational wave packet (coherent superposition of rotational states) in the molecule. By 
treating the linear molecule as a rigid rotor j^, [3] , the rotational motion of the molecule 
with initial state jm{0, (f),t = — oo) = \ JM) evolves in the laser field following the time- 
dependent Schrodinger equation 

^ d^jMM,t) ^ ^^j, _ Ei^^^^^ ^ ^ e)]^jM{e, <p, t). (47) 

Here E{t) is the laser electric field, B is the rotational constant, a\\ and a± are the anisotropic 
polarizabilities in parallel and perpendicular directions with respect to the molecular axis, 
respectively. These molecular properties for CO2, O2 and N2 are given in Table. I. The above 
equation is then solved for each initial rotational state | JM) using the split-operator method 
[see Eq. ([8])]. We assume the Boltzmann distribution of the rotational levels at the initial 
time. With this assumption, the time-dependent alignment distribution can be obtained as 

p{9,t)=Y,u;jM\'^jM{e,M\', (48) 

JM 

where oojm is the weight according to the Boltzmann distribution. Note that one needs to 
take proper account for the nuclear statistics and symmetry of the total electronic wave- 
function. For example, in case of O2 with total electron wavefunction in only odd-J's 



9 

are allowed (see for example 73|), whereas CO2 with total electron wavefunction in 



g 

has only even- J's. The angular distribution or alignment does not depend on the azimuthal 
angle in the frame attached to the pump laser field. The two equations above allow the 
determination of the time- dependent alignment distribution of the molecules in the laser 
field, as well as the rotational revivals after the laser has been turned off. The aligning laser 
is assumed to be weak enough so the molecules remain in the ground state and no ionization 
occurs. 

Once the angular distribution is obtained, the (complex) induced dipole for emission of 
photon energy of u can be calculated by adding coherently the weighted contribution from 
different alignments by 

D(uj, t) = 27T r D{uj, e)p{e, t) sin OdO, (49) 

JO 
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TABLE I: Molecular properties for CO2, O2 and N2. B is rotational constant, ay and a± are 
parallel and perpendicular polar izability, respectively. The data are taken from 



74 



Molecule 


B (cm-i) 


«ll (A3) 


a± (A3) 


CO2 


0.39 


4.05 


1.95 


O2 


1.4377 


2.35 


1.21 


N2 


1.989 


2.38 


1.45 



if the pump and probe laser polarizations are parallel. Here we assume that rotational motion 
during the femtosecond probe pulse is negligible, which should be valid for the molecules 
under consideration with typical rotational period of few picoseconds. 

If the polarizations of the pump and probe lasers are not the same, the theoretical treat- 
ment is rather cumbersome as the cylindrical symmetry is lost. This case has been discussed 
by Lein et al j?^. Here we briefly describe the main results. Assume that the pump and 
probe laser pulses propagate collinearly and a is the angle between the two polarization 
directions. Let 9 {9') and (0') be the polar and azimuthal angles of the molecular axis in 
the frame attached to the pump (probe) field. These angles are related by 

cos 9 = cos 9' cos a + sin 9' sin a cos 0'. (50) 
The alignment distribution in the "probe" frame is 

pia,9',<P',t) = p{9ia,9',<P'),t). (51) 

For the emitted HHG with polarization parallel to that of the probe laser, the induced dipole 
can then be obtained from 

Dn (lu, a, t) = r r Dn {lu, 9')p{a, 9', 0', t) sin 9'd9'd(j)'. (52) 
Jo Jo 

And for the perpendicular component 



D^{uj,a,t) = r f " D^{uj,9',(j)')p{a,9',(f)',t)sm9'd9'-'^' 
Jo Jo 

= r D±{uj, 9', (j)' = 0)p{a, 9', (f)',t) sin 9' cos (j)'d9'd(f)'. (53) 
Jo Jo 
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III. QUANTITATIVE RESCATERING THEORY 



In this section we provide formulation of the quantitative rescatering theory (QRS) and 
theoretical evidence in supporting its validity and improvements over the very popular SFA 
(or the Lewenstein model). Particular attention is given to the photo-recombination dif- 
ferential cross sections (DCS) and phases, which are used to illustrate the nature of the 
improvements and will also be used in the next section to simulate data for comparison with 
experiments. 



A. Description of the QRS 

Within the QRS as applied to HHG process, induced dipole D{uj,9) and its phase {p{uj,9) 
for a molecule aligned with an angle 6 with respect to the laser polarization can be written 

as 

D{uj,e) = w{E,e)d{uj,e), (54) 

or more explicitly 

|D(cu,^)|e''^("'^) = \ W{E,9)\e'''\d{iu,e)\e'^^'''^\ (55) 

where d{uj, 9) and 5{uj,9) are the "exact" transition dipole and its phase defined in sections 
11(C) and (D). The quantity |W^(-E,6')p describes the flux of the returning electrons, which 
we will call a "wave packet", with r]{E,9) being its phase. Electron energy E is related 
to the emitted photon energy uj hj E = uj — Ip, with Ip being the ionization potential of 
the target. Clearly the HHG signal S{uj,9) ~ \D{uj,9)\'^ and W{E,9) depend on the laser 
properties. On the other hand, d{u!, 9) is the property of the target only. Equation (154|) can 
be also seen as the definition of the wave packet, assuming the induced dipole and transition 
dipole are known. Since the returning wave packet is an important concept in the QRS 
theory, let us write it down explicitly 



The validity of Eq. (jSlD on the level of amplitudes has been shown in Morishita et al. 38 1 



using HHG spectra calculated by solving the TDSE for atoms. Indications for the validity 



of this factorization have also been shown for rare gas atoms by Levesque et al. 77|] and 
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for N2 and O2 molecules [3] , where the HHG spectra were calculated using the SFA model 
with the continuum electron being treated in the plane wave approximation. 

In the tunneling regime, most electrons will be driven along the laser polarization direc- 
tion. In fact, semi-classical treatment shows that most contribution to the HHG process is 
coming from electron released to the continuum and returning to the target ion along the 
laser polarization direction {22]. In the notations adopted in Sec. 11(C) and (D), with n and 
n' being the directions of the driving laser and HHG polarizations, respectively, and k the 
electron momentum, we have 

{fc II n II n' for -Diifu;, 9) 
liv , ; ^g^^ 
k\\n±n' ioT D±{uj,e). 

In this paper, we limit ourselves to the parallel component of HHG, and therefore the 
subscripts are omitted in the notations. 

The usefulness of Eq. is two-fold. First, the factorization allows us to separate in 
HHG process the effect of the laser field (on the returning electron wave packets) and the 
influence of target structure (transition dipole and its phase). From practical point of view, 
this implies that one can carry out calculations for each factor separately. It is particularly 
important to be able to use the vast knowledge of the molecular photoionization/photo- 
recombination processes, which have been accumulated over the last few decades. As for 
the wave packets, it will be shown that although the HHG spectra from the SFA are not 
quite accurate, the extracted wave packets are reasonably good. This offers a simple and 
efficient way to calculate the wave packets. Second, the wave packets can also be shown to 
be largely independent of the targets. By that we mean that the shape of the wave packet 
as a function of electron energy depends only on laser parameters, for targets with similar 
ionization potentials. In other words, the wave packets can be written as 

\W{E,9)\' = N{9,...)x\wm', (58) 

where N{6,...) is the ionization probability for the emission along the laser polarization 
direction, which depends on the alignment angle, symmetry of the HOMO, and other pa- 
rameters. The overall factor N{6, ...) does not change the shape of W{E). This offers another 
way to obtain the wave packets by using a reference atom, for which numerical solution of 
the TDSE is relatively simpler than that of a molecule under consideration. Regarding the 
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factor N(6, ...), it is well understood for many molecular systems based on the MO-ADK 



theory 
systems 



67|and MO-SFA [see Sec. 11(E)], and the numerical solution of the TDSE for simple 



26 



27 



79 



80[. Note that, the numerical solution of the TDSE for ionization rate 



is much less computational demanding than for calculating HHG spectra. 

From the above general discussion, let us now be more specific about the two ways of 
obtaining the wave packets and HHG spectra within the QRS model. For a given target, we 
can obtain the wave packets from the SFA. If the laser pulse is given, one first uses the SFA 
to calculate D^^^(t, 0) and its Fourier transform D^^^{u, 6) by using Eq. (11). Taking into 
account that the PWA is used in the SFA, we have 

D^^^{uj, 6) 



(59) 



It will be shown that this wave packet agrees reasonably well with the "exact" wave packet 
W{E, 9) obtained from solving the TDSE. Once the wave packet is obtained, the induced 
dipole can be calculated by 



d{uJ,9) ^SFA, 



u,9). 



(60) 



dP^^{uj,9) 

Since the photo-recombination cross section is proportional to modulus square of the tran- 
sition dipole [see Eqs. (14), (15), (25), and (33)], the HHG yield can be written as 



d{Lj,9) 



dP^\uj, 9) 
aiuj, 9^ 



a 



PWA, 



u;,9) 



S''^^{uj,9) 
S'''\uj,9), 



(61) 



where a and a^^^ are the short-hand notations for the "exact" and PWA differential photo- 
recombination cross sections, respectively. For simplicity, we have used QRSl to denote 
this version of the QRS model. Since the wave packets are obtained from the SFA, this 
method was called the scattering-wave based strong-field approximation (SW-SFA) in our 



previous papers 



36, 



371 ]. This method has a very simple interpretation: the QRS corrects 



the inaccuracies in the HHG yield from SFA by a simple scaling factor, equal to the ratio 
of the exact and approximate photo-recombination cross sections. It also adds the exact 
transition dipole phase to the harmonic phase. The success of this version of the QRS lies 
in the fact that the SFA describes the electron wave packet quite accurately. Recall in the 



21 



three-step model, in the second step the electron "roams" well outside the target ion most 
of the time before being driven back to recoUide with the ion. During this excursion, the 
electron motion is governed mostly by the laser field, which is well described by the SFA. 
Thus, in this version only the transition dipole moment is corrected. 

The second method of obtaining the wave packet for the QRS is to use a reference atom 
with a similar ionization potential. For the reference atom, we can perform the TDSE 
calculation. Using the idea of the QRS, one can obtain the wave packet from 

^''^^^E) = (62) 

The power of this method stems from the fact that effect of the target potential on the 
wave packet is included to some extent in the second step of the three-step model (see the 
paragraph above) , when electron is quite far from the ion core and sees mostly the long-range 
Coulomb tail. This method has the advantage in improving the accuracy of the phase of 
the HHG induced dipole, but it is much more time consuming. Combining with Eq. (67), 
the wave packet for the molecular system of interest can then be written as 

1/2 



where and N'^^^ are the ionization probabihty for electron emission along the laser 

polarization direction from the molecule and reference atom, respectively. A77 is introduced 
to account for the phase difference between the two wave packets. This phase difference will 
be shown to be nearly independent of energy. The induced dipole and HHG spectra can 
then be written as 

^ (N{e)V'' d{u,e) 



and 



TTT TTP^^''"'D-f{ul (64) 

SQ^''{u,e) = ^^^^S-fico). (65) 
^ ^ N^'^f (7''^/ (a;) 

This version of the QRS, called QRS2 above, also has important implications. It reflects 
the fact that in the tunneling regime the returning electron wave packet has nearly identical 
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shape (momentum distribution) for all targets, except for an overall factor accounting for the 
differences in ionization rates. In practice for the reference atom we use a scaled hydrogen 
with the effective nuclear charge chosen such that it has the same Is binding energy as the 
molecule under consideration. Experimentally, one can replace the scaled atomic hydrogen 
with an atomic target of comparable ionization potential. The wave packet obtained from 
the reference atom also has the advantage that it avoids the spurious singularity often seen 
in the wave packet obtained from the SFA. Such singularity occurs since the transition 
dipole calculated from PWA usually goes to zero at some photon energy, see Eq. (!60l) . Both 
versions of the QRS have been used. In general, the one using the reference atom is more 
accurate but takes much longer time since the wave packets are obtained from the TDSE 



36 



37|, 



38 



8l| for other 



for atoms. Interested reader is also referred to our recent papers 
evidences and applications of the QRS to HHG processes. We note that the factorization 
can be approximately derived analytically based on the SFA |t82| and within the zero-range 
potential approach [83 1. 



B. Atomic photo-recombination cross sections and phases 



In this subsection we will consider an example of photo-recombination of Ar(3po), which 
will be used in the QRS calculation in the next subsection. We will use a single-active 
electron model with a model potential, as detailed in Sec. 11(A) and (C). We have found 
that the position of the Cooper minimum is quite sensitive to the form of the model potential. 



In our earlier paper 



36 



38| we used the model potential of Tong and Lin [50|, which shows 



a Cooper minimum near 42 eV. To have more realistic simulations for HHG spectra in this 

n 

paper we use a model potential, suggested by Muller [84], which has been shown to be able 
to reproduce the ATI spectra comparable with experiments [85] • This model potential has 
also been used quite recently to simulate HHG experiments by Minemoto et al 86|] and 



Worner et al 



871 ]. In Fig. 1(a) and 1(b) we compare the differential photo-recombination 



cross sections and dipole phases obtained from Eqs. fl25l) and (!33l) with these two model 
potentials. We also plot here the results from the PWA, with the ground state wavefunction 
from Muller potential, which are almost identical with the PWA results from Tong and 
Lin potential (not shown). Clearly, in the energy range shown in the figure the PWA 
result deviates significantly from the two more accurate results. Furthermore, the Cooper 
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FIG. 1: (Color online) Comparison of the differential photo-recombination cross sections (a) and 
dipole phases (b) from different models for Ar. 

minimum obtained from the potential by Muller occurs near 50 eV, in a better agreement 
with experiments, although the minimum is somewhat more shallow than that from Tong 
and Lin model potential. This is not surprising since the scattering waves are known to be 
more sensitive to the details of the potential than the bound states, which are localized near 
the target core. The dipole phases obtained from both potentials show dramatic jumps of 
about 2 radians near the Cooper minima. Note that the dipole phase from the PWA shows 
a phase jump by vr at the "Cooper minimum" near 21 eV. The failure of PWA for describing 
photo-recombination cross section and phase at low ener gies is well-known. Other examples 
can be found in Ref. 36| for rare gas atoms and in Ref. [37!] for H^. 

The simplicity of the model potential approach allows us to establish the validity of the 
QRS for atoms, by using the same potential in the TDSE and in the calculation of the photo- 
recombination transition dipole (see next subsection). In order to compare with experiments, 
one can go beyond the model potential approach by using different schemes to account for 
many-electron effect, i.e., final state and initial state correlation in photoionization (see, for 



example. 



56 



88 



891]). It is much harder to do so within the TDSE approach. 
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FIG. 2: (Color online) (a) Comparison of the returning electron "wave packets" extracted from 
numerical solutions of the TDSE for Ar and scaled hydrogen and from the SFA for Ar. (b) 
Comparison of the HHG yields obtained from numerical solutions of the TDSE, QRS and SFA for 
Ar. Data have been shifted vertically for clarity. 8-cycle laser pulse with peak intensity of 2.5 x 10^'^ 
W/cm^, 800- nm wavelength is used. 

C. QRS for atomic targets: example of Argon 

By using the photo-recombination cross section and dipole phase shown in the previous 
subsection, we are now ready to discuss the results from the QRS and assess its validity 
on an example of HHG from Ar(3po). We emphasize that in this paper we use a model 



potential proposed by MuUer 8^ in order to get the Cooper minimum position comparable 
with experiments. 

In Fig. 2(a) we show the comparison of the returning electron wave packets from both 
versions of the QRS with the "exact" wave packet, extracted from the solution of the TDSE 
for Ar by using Eq. (65). In the QRS2 the effective charge Zejf = 1.0763 is chosen such that 
the ionization potential from Is state is 15.76 eV, the same as for Ar(3po). We use a 800-nm 
wavelength laser pulse with 8-cycle duration (8 fs FWHM) and peak intensity of 2.5 x 10^^ 
W/cm^. One can see very good agreement between QRS2 and the exact result over a very 
broad range of energy when the wave packet is extracted from the scaled atomic hydrogen. 
The result from the QRSl (shifted vertically for clarity) with the wave packet from the SFA 
is also in reasonable good agreement with the exact one. 

Fig. 2(b) shows comparison of the HHG yields from the TDSE, QRSl, QRS2 and the 
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FIG. 3: (Color online) (a) Comparison of extracted harmonic phase difference Ac/? between Ar and 
scaled hydrogen as a function of emitted photon energy with the photo-recombination transition 
dipole phase difference A6. The results are obtained from numerical solutions of the TDSE. (b) 
Same as for (a), but SFA for Ar is used instead of a scaled hydrogen. Laser parameters are the 
same as for Fig. 2. 

SFA. The data from QRSl and SFA have also been shifted vertically for clarity. Clearly, the 
QRS results agree quite well with the TDSE, whereas the SFA result deviates strongly in 
the lower plateau. The signature of the Cooper minimum near 50 eV in the HHG spectra 
is quite visible. Note that the minimum has shifted as compared to the results reported 



in Ref. 36|, which used a model potential suggested by Tong and Lin [5(1]. The results in 
Fig. 2(b) clearly demonstrate the good improvement of the QRS over the SFA in achieving 
better agreement with the TDSE results. 

Next we examine harmonic phase, or induced dipole phase ip [see Eq. fl55|) ]. Since har- 
monic phase changes quickly from one order to the next (except in the harmonic cut-off), it 
is instructive to compare phases from different targets. For that purpose we calculate the 
harmonic phase difference Aip = ip — ip^'^-^ of Ar from its reference scaled hydrogen. The 
result, shifted by 1.9 radians, is shown as the solid black curve in Fig. 3(a), at an energy 
grid with step-size of 0.1 x ujq, where ujq is the photon energy of the driving laser (1.55 
eV). A striking feature is that this curve agrees very well with the transition dipole phase 
difference A6 = 6 — 6^'^^ , shown as the red line. Since the transition dipole phase S^'^'^ from 
the reference scaled hydrogen is very small (about 0.2 radians) in this range of energy, A5 
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FIG. 4: (Color online) The electron wave packets, extracted from the HHG spectra for Ne and Ar 
under the same laser field. The inset shows the wave packets in the linear scale with the Ne data 
shifted horizontally by -1.2 eV. The HHG are obtained within the Lewenstein model with 1064 nm 
laser, peak intensity of 2 x 10^^ W/cm^, duration of 50 fs. 

looks quite close to 6 as well [see Fig. 1(b)]. Note the phase jump near the Cooper minimum 
at 50 eV is also well reproduced. This clearly demonstrates the validity of the QRS model in 
Eqs. ( 155|) and ( IMl) and with respect to the phase. The shift by 1.9 radians can be attributed 
to the phase difference of the two wave packets Arj, which is nearly energy independent. 

Above result is based on the numerical solution of the TDSE for Ar and scaled hydrogen. 
Similarly one can compare the harmonic phase difference between the TDSE and SFA results 
= Lf-Lp^^"^ with that of the transition dipole phase AS = Note that S^^"^ = 

[see Fig. 2(b)]. Clearly, if the SFA and the PWA were exact, one would have Aip = A6 = 0. 
The results, presented in Fig. 3(b), show that the phase missed in the SFA is quite close to 
the transition dipole phase. The QRSl therefore corrects the inaccuracy in the SFA phase 
by adding the phase from the transition dipole [see Eq. (!60l) ]. By comparing Figs. 3(a) 
with 3(b), we conclude that the phase from QRS2 is more accurate than QRSl. This is 
expected since the QRS2 includes partially the effect of electron-core interaction during the 
propagation in the continuum, as explained in III (A). 

Finally, we remark that the use of a reference atom with nearly identical ionization 
potential is not necessary. We will show here this requirement can be relaxed. As the 
wave packets obtained within the SFA agree reasonably well with the TDSE results, for our 
purpose we only use the SFA in the subsequent analysis. In Fig. 4, we show the comparison 



27 



of the wave packet from Ne {Ip = 21.56 eV) and Ar (Jp = 15.76 eV) in the same 1064 nm 
laser pulse {uj = 1.166 eV) with peak intensity of 2 x lO^'' W/cm^, duration (FWHM) of 50 
fs. Note that the wave packets are now plotted as functions of electron energy, instead of 
the photon energy as in Fig. 2(b). Clearly, the two wave packets lie nicely within a common 
envelope. This can be seen even more clearly in the inset, where the data are plotted in linear 
scale with the Ne data shifted horizontally by —1.2 eV. We note that the small details below 
33 eV also agree well. This conclusion is also confirmed by our calculations with different 
laser parameters and with other atoms. This supports that the independence of the wave 
packet on the target structure can also be extended to systems with different ionization 
potentials. This fact can be useful in comparing experimental data for different targets. 
Note that the small shift of -1.2 eV is caused by the fact that the difference in the ionization 
potentials (5.8 eV) is incommensurate with the energy of two fundamental photons (2.3 eV). 

D. Molecular photoionization cross sections and phases 

Photoionization cross section and transition dipole phase from a linear molecule depend 
on both photon energy and angle between molecular axis and laser polarization direction. 
Since most interesting features of HHG from molecules are from molecules that have been 
pre-aligned by a pump laser pulse, it is more appropriate to present them as functions of 
fixed alignment angle for fixed energies. As noted before, one can use either photoionization 
or photo-recombination cross sections and phases in the QRS. Here we choose to use pho- 
toionization as it is in general more widely available theoretically and experimentally. In 
this paper we are only interested in the HHG with polarization parallel to that of the driving 
laser. That limits the differential cross sections to that case of || n || n' [see Eq. f l57|) ]. 

The photoionization DCS and phase from O2, CO2, and N2 are presented in Fig. 5(a-f) 
for some fixed energies as functions of angle between molecular axis and laser polarization 
direction. For convenience we express energy in units of photon energy of 800-nm laser. Let 
us first examine O2, shown in Fig. 5(a) and (b). For a fixed energy the cross section vanishes 
at = and tt/2 due to the iig symmetry of the HOMO and the dipole selection rule for the 
final state and has a peak near 45°. As energy increases the peak slightly shifts to a larger 
angle and the cross section monotonically decreases. The phases behave quite smoothly and 
change only within about 1 radian for all energies considered. For CO2, the cross section 
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FIG. 5: (Color online) Diffrential photoionization cross sections and transition phases for O2, CO2, 
and N2, as functions of the alignment angle. The photon energies are expressed in units of harmonic 
orders for 800-nm laser. 
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also vanishes at ^ = and 7r/2 due to the vr^ symmetry of the HOMO. However in contrast 
to O2, for a similar range of energy the DCS of CO2 shows a double- hump structure with the 
minimum shifting to a larger angle as energy increases. These features can also be seen in 
the phase of the transition dipoles [Fig. 5(d)], where the phase jump by almost tt is observed 
near the minima in the cross section. The two "humps" also behave quite differently. Within 
the range of energy from H21 to H39, the small hump at small angles increases and the big 
hump at large angles decreases with increasing energy. We note that the phase jump is 
smaller than vr as the cross section does not go to zero at the minimum between the two 
humps. In this paper we will limit ourselves to application of the QRS to O2 and CO2. 
Nevertheless it is worthwhile to point out that the case of N2 is even more complicated 
than CO2, with almost no regular behavior found in the energy range presented here, see 
Fig. 5(e) and (f). Note that the sharp variation of cross section from H17 to H25 is due to 
the presence of the well-known shape resonance in N2. The shape resonance is accompanied 
by a rapid phase change in the same energy region. 

We have seen that these three molecular systems behave totally differently in photoion- 
ization. We show in Sec. IV that such differences can be seen from the HHG spectra. 



E. Wave packets from molecular targets: example of O2 

In subsection III(C) we have shown evidence to validate the QRS for atomic targets. For 



molecular targets some benchmark TDSE results are available only for H2 [13|, |26|, |27|]. In 
fact successful application of the QRS and detailed comparisons with the TDSE results for 
this system has been reported in js?]. In this subsection for completeness we show here 
an example of the wave packets from O2 and Xe, which have nearly identical ionization 
potentials (12.03 eV for O2 and 12.13 eV for Xe). In Fig. 6(a) we compare the HHG spectra 
from O2 aligned at 10° and 70°, with that from Xe. These results were obtained from the 
SFA, with 1600 nm laser of 1 x lO^'' W/cm^ intensity and 20 fs duration (FWHM). The long 
wavelength is used here in order to compare the wave packets in the extended HHG plateau. 
Clearly, the HHG spectra look quite different. However, the wave packets, extracted by 
using Eqs. (!59l) and ( |62l) . shown in Fig. 6(b) are almost identical. For clarity we have shifted 
the curves vertically in both Figs. 6(a) and (b). The QRS results for the HHG are shown 
in Fig. 6(c). The changes in the slope of the QRS yields as compared to that of the SFA 
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in Fig. 6(a) can be easily seen, which reflects the differences of the exact and the PWA 
photoionization DCS. 

This example clearly shows that in order to calculate the HHG spectra for a flxed align- 
ment within the QRS one can use the wave packet from any other alignment angle or from a 
reference atom, with an overall factor accounting for the differences in tunneling ionization 
rates [see Eq. (l58l) ]. This feature enables us to avoid the possible difficulties associated with 
direct extraction of the wave packet at the energies and alignments where the transition 
dipole in the PWA vanishes. For example, in case of N2, it is more convenient to use the 
wave packet extracted for alignment angle of 90°, as the transition dipole in the PWA does 
not vanish. Furthermore, if a reference atom is used, one can also extract the wave packet 
from the solution of the TDSE, as has been done for js^]- 



IV. COMPARISONS WITH EXPERIMENTS 



In this section we will compare the QRS results with the recent HHG measurement from 
partially aligned molecules. Below we will look at two examples of molecular O2 (with 
internuclear distance at equilibrium R = 2.28 a.u.) and CO2 {R = 4.38 a.u. between the 
two oxygen centers), which can also be seen as "elongated" O2. HHG from aligned molecules 
are typically measured experimentally with the pump-probe scheme. In this scheme, non- 
adiabatic alignment of molecules is achieved by exposing them to a short and relatively 
weak laser pulse (the pump) to create a rotational wave packet. This wave packet rephases 
after the pulse is over and the molecules are strongly aligned and anti-aligned periodically 
at intervals separated by their fundamental rotational period 9|. To observe the alignment 
dependence of HHG, a second short laser pulse (the probe) is then used to produce HHG 
at different short intervals when the molecules undergo rapid change in their alignment. A 
slightly different setup is done by changing the relative angle between the pump and probe 
polarizations, but measure the HHG at fixed time delays, typically at maximal alignment 
(or anti-alignment) near half-revival. Theoretically, induced dipoles from the QRS need to 
be calculated first for a fixed molecular alignment. The results will then be convoluted with 
the molecular distribution following the theory presented in Sec. 11(F). 
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FIG. 6: (Color online) (a) Comparison of the HHG spectra obtained from the SFA for O2 aligned 
at 10° and 70°, and Xe under laser pulse with wavelength of 1600 nm, intensity of 1 x 10^^ W/cm^ 
and duration of 20 fs. For clarity, we have shifted Xe data vertically, (b) Comparison of the 
electron returning wave packets extracted from (a). Wave packet for O2 at 10° is renormalized to 
that of 70° and Xe and plotted twice, (c) The HHG spectra from O2 aligned at 10° and 70° from 
the QRS. 
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A. HHG from aligned O2 molecules 



In Sec. ni(D) we showed harmonic spectra generated with 1600-nm laser. This was used 
in order to compare the wave packets in an extended plateau. To compare with available ex- 
periments calculations in this subsection are performed with 800-nm laser. First we show in 
Fig. 7(a) the HHG yields for some selected harmonics for fixed molecular axis. Calculations 
are done with the QRS, for a 30-fs pulse with a peak intensity of 2 x 10^'* W/cm^. The yields 
are maximal if molecules are aligned at about 45° with respect to the probe polarization, 
and vanish at 0° and 90° due to vr-symmetry of the HOMO. This result resembles closely 
the behavior of the photoionization DCS's shown in Fig. 5(a), but with the relative magni- 
tudes changed, reflecting the influence of the returning wave packets. The yields convoluted 
with the partial distribution with maximally aligned ensemble at half-revival are shown in 
Fig. 7(b), as functions of relative angle between pump and probe polarizations. Interest- 
ingly, the yields are now peaked near 0°, that is, when the pump ar id p robe polarizations are 
parallel. This result is consistent with the data by Mairesse et al 901] (see their Fig. 2) and 



with the earlier measurements by Miyazaki et al [9l| . Note that the yield is quite insensitive 
to the alignment angle. However, the contrast can be increased with a better alignment. 
Our simulations for alignment are carried out using the pump laser intensity of 5 x 10^^ 
W/cm^ and 30 fs duration, with a rotational temperature of 30 K, the same as used in the 
experiment. We note that the convoluted yields from the earlier SFA results by Zhou el al 
291] show peaks near 45° with similar laser parameters. Within the QRS we found that only 
with higher degrees of alignment the peaks shift to about 45°. This can be achieved, for 
example, with the same pump laser intensity, but with a longer duration. 

Next we compare in Fig. 8 harmonic yields from H17-H35 as a function of delay time near 
quarter- and half-revivals from the QRS and the SFA calculations with the experiments by 



10| and 



Itatani et al 1^. Calculations were done with laser parameters taken from Ref. 
the yields have been normalized to that of the isotropic distribution. Clearly both theories 
reproduce quite well the behavior of the experimental curve, which follows < sin^2^ >, 
shown in the upper panel. A closer look reveals that the QRS result agrees better with 
the experiment and the SFA tends to overestimate the yield near maximal alignment and 
underestimates it near maximal anti-alignment. In fact, similar quantitative discrepancies 
between the SFA and experiments can be seen in an earlier work by Madsen et al 92| for both 
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Alignment angle (degrees) Angle between pump & probe polarizations (degrees) 



FIG. 7: (Color online) Alignment dependence of harmonic yields from O2 harmonics H21, H25,H29, 
and H33 for fixed molecular axis (a) and convoluted with molecular alignment distribution (b). 




FIG. 8: (Color online) Comparison of O2 time evolution of (normalized) harmonic yields H17-H35 
from the QRS, SFA and experiment, near quarter and half-revivals. The experimetal data and 
laser parameters are taken from Itatani et al [l^. < sin^ 2^ > is shown in the upper panel for 
reference. 

O2 and N2. Note that the discrepancies between the simulations and the experiment along 
the delay time axis could be due to inaccuracy in the determination of laser intensities and 
temperature used in the experiment. We have seen so far only quantitative improvements 
of QRS over the SFA for O2. The situation is completely different for CO2, presented in the 
next subsection, where the two theories predict qualitatively different results. 
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B. HHG from aligned CO2 molecules 



First we show in Fig. 9(a) theoretical HHG yields for selected harmonics from H21 up to 
H39, obtained from the QRS for fixed angles between the molecular axis and the polarization 
of the probe pulse. The yields resemble closely the photoionizations DCS of Fig. 5(c), with 
the double-hump structures seen quite clearly for the higher harmonics above H31. For the 
lower harmonics, the smaller humps nearly disappear. To compare with experiment, average 
over the molecular alignment needs to be performed. We take the alignment distribution 
for a pump-probe delay time at the maximum alignment near half-revival. The convoluted 
yields are presented in Fig. 9(b), as functions of the angle between the pump and probe 
polarizations. Clearly, due to the averaging over the molecular alignment distributions, the 
angular dependence of HHG is smoother as compared to fixed alignment data in Fig. 9(a). 



These results are consistent with recent experiments |18l . l90l . |93| . which show enhanced 
yields for large alignment angles and minima near 30° for harmonics above H31. We note 
that our results also resemble the data for the induced dipole retrieved from mixed gases 
experiments by Wagner et al. (see their Fig. 4). To have a more complete picture of the 
HHG yields, we show in Fig. 9(c) a false-color plot of HHG yield as a function of harmonic 
orders and the angle between the pump and probe lasers polarization directions, at the delay 
time corresponding to the maximum alignment near the half-revival. Clearly the yield has 
quite pronounced peak at large alignment angle. However, the most pronounced feature 
in Fig. 9(c) is a minimum at small angles, which goes to larger angles as with increasing 
harmonic orders. Within the QRS the origin of this minimum can be directly traced back 
as due to the minimum in the photo-recombination DCS's. Our results are comparable with 
the recent measurement by Mairesse et al 90| (see their Fig. 2), which were carried out 



with a slightly different laser parameters. Interestingly, our results resemble closely with 



the theoretical calculations by Smirnova et al [9J], which include contributions from the 
HOMO, HOMO-1 and HOMO-2 orbitals [see their Fig. 1(b), plotted in a x-y plot]. In our 
calculations, only the HOMO is included. We note that the position of the minima can be 
slightly shifted depending on the actual experimental laser setup. The above calculations for 
HHG were done with a 25-fs probe laser pulse with intensity of 2.5 x 10^^ W/ cm^. Alignment 
distribution was obtained following the method is Sec. 11(F), with a 120-fs pump laser pulse 
with intensity of 0.55 x 10^^ W/cm^. Rotational temperature is taken to be 105 K. These 
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parameters are taken from the recent experimental setup by Zhou et al [l^ . 

Comparison of theoretical HHG amplitudes for a few angles between pump and probe 
polarizations from 0° to 25° is shown in Fig. 10(a) together with the experimental data in 
Fig. 10(b), taken from Mairesse et al 90|]. The experimental data have been renormalized to 
a smoothed experimental wave packet from Ar under the same laser field. This wave packet is 
a decreasing function of harmonic order, therefore the renormalized signals in Fig. 10(b) are 
enhanced at higher orders compared to lower orders. No such renormalization was done for 
the theoretical data. Nevertheless, reasonably good agreement with the experimental data 
can be seen, including the shift of the minimum position to higher harmonics as alignment 
angle increases. The theoretical data seem to have more structures than the experimental 
ones. That might be due to the fact that that we use a single intensity, single molecule 
simulation. It is known that macroscopic propagation in general tends to smooth out the 
HHG spectra. 



As found in Zhou et al [16|], the most dramatic feature can be seen near 3/4- revival, 
where the molecule can be most strongly aligned. Our QRS results are shown in Fig. 11 for 
the same harmonics that have been analyzed in IGj (see their Fig. 2). (We note here that 
according to the authors of [l^, the harmonic orders should be properly shifted down by 
two harmonic orders, as compared to the ones given in their original paper). Theoretical 
results were earned out with >ase. pata,nete. taken .on, As can be seen from our 

results for the lower harmonics H21 and H25, the yields follow the inverse of the alignment 
parameter (cos^O), shown in Fig. 11(a). However, for the higher harmonics, an additional 
peak appears right at the delay time corresponding to the maximum alignment. The peak 
starts to appear near H31 and gets more pronounced with increasing harmonic orders. This 
behavior is in quantitatively good agreement with the measurements by Zhou et al. 
In Zhou et al, the experimental data were fitted to the two-center interference model, by 
using a three-parameter least-square fitting procedure. Within the QRS, no fitting is needed. 
Furthermore, one can trace back the origin of the time-delay behavior based on the two-hump 
structure of the photoionization DCS. Indeed, Fig. 5(b) shows that for the lower harmonics 
H21 and H25 the cross sections at large angles are much larger than at small angles so the 
HHG yield is inverted with respect to (cos^O), for which small angles dominate. This fact 
has been known before [sO^. For the higher orders. Fig. 5(b) indicates that the smaller hump 
at small angles become increasingly important. Qualitatively, this explains why the HHG 
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FIG. 9: (Color online) HHG yield as a function of fixed angle between molecular axis and laser 
polarization (a) and as a function of angle between pump and probe laser polarization directions 
(b). False-color plot of HHG yield for a harmonic range between H20 and H42 (c). In (c) the probe 
laser is horizontally polarized. Alignment distribution is chosen at the maximum alignment near 
half-revival. The laser intensity and duration are of 0.55 X 10^'' W/cm^, and 120 fs for the pump 
pulse, and 2.5 x 10^^ W/cm^ and 25 fs for the probe pulse. Rotational temperature is taken to be 
105 K. 
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amplitudes for CO2 for a few angles from 0° to 25°. The 



FIG. 10: (Color online) Harmonic 
experimental data by Mairesse et al 90(] are also shown for comparison (b). The experimental data 
have been renormalized to a smoothed experimental wave packet from Ar under the same laser 
field. 

yields for H31 and up show a pronounced peak for the parallel alignment. In fact, it is 
even more simpler to explain this behavior by looking directly at the HHG yields shown in 
Fig. 9(a). One can immediately notice that the peak at small angles near 25° starts to show 
up more clearly only near H31. We note that apart from the additional peaks, which are 
clearly visible at 3/4- revival, the general behavior of all harmonics from H19 up to cutoff 
near H43 shows inverted modulation with respect to (cos^^). Therefore, in contrast to the 
interpretation of Kanai et al 1^ , inverted modulation is not an unambiguous indication of 
the interference minimum, especially for determining the precise position of the minimum. 
Since the laser intensity in the experiment is quite high, the depletion of the ground state 



needs to be accounted for 
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951]. Therefore discussion about the ionization rate is in 



order. Alignment dependence of the ionization rate for CO2 is still a subject of debate. In our 



previous attempts 



30 



3l| , the MO-ADK theory 67|] has been used to calculate the ionization 



rate. The MO-ADK theory predicts the peak in ionization rate near 30°, in agreement with 
the results deduced from measured double ionization 96|. The recent experiment by the 
NRG group [71], however, show very narrow peak near 45°. Interestingly, the MO-SFA 
predicts a peak near 40°. However the MO-SFA theory is known to underestimate the 
ionization rate. In order to correct the MO-SFA rate, we renormalize the MO-SFA rate 
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FIG. 11: (Color online) Normalized HHG yield (vs isotropically distributed molecules) from CO2 
for selected harmonics as function of pump-probe delay time near 3/4-revival. The experimental 
data are taken from Zhou et al jl6]. Theoretical results were obtained with laser parameters taken 
from [m]. The alignment parameter (cos^ 0) is also plotted for reference in the top panel. 

to that of the MO-ADK rate at laser intensity of 10^^ W/cm^, which gives a factor of 10. 
Note that the same correction factor has been found for the SFA ionization from Kr, which 
has almost the same ionization potential as for CO2. With the corrected MO-SFA rate, we 
found that our simulations give a better quantitative agreement with the JILA data, than 
with the MO-ADK rate. On the other hand, calculation with the uncorrected rate would 
lead to much more pronounced peaks right at the maximum alignment near 32.6 ps for H31 
and higher harmonics. 

We next compare the predictions of the QRS for harmonic phase with experiments. The 
phase information can be of great importance. In particular it has been suggested as a more 
accurate way to determine positions of interference minima in HHG spectra. Experimentally 



the harmonics phase can be extracted from measurements of HHG using mixed gases 18l.ll9| 
or interferometry method [igI. In Fig. 12(a), we show the recent experimental data of Boutu 



et al. [18| where the harmonics phases (relative to that from Kr) are obtained for the parallel 
aligned and perpendicularly aligned ensembles, shown by black and red symbols, respectively. 
For the latter, the phase does not change much within H17 to H31. For the parallel aligned 
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FIG. 12: (Color online) (a) Harmonic phase (relative to the phase from Kr) for parallel aligned 
and perpendicular aligned ensembles, under two sets of parameters (solid and dashed lines) that 
lead to two different degrees of alignment distributions. Experimental data from Boutu el al 
are shown as symbols. Calculations are carried out with the same probe laser parameters as in the 
experiments, (b) Phase of the transition dipole for alignment angles from 20° to 60°. 

molecules, the phase jump from H17 to H31 was reported to be 2.0 ± 0.6 radians. Our 
simulations are shown for two different ensembles, with alignment distributions confined 
in a cone angle of 25° (solid lines) and 35° (dashed lines) at half maximum. The less 
aligned ensemble was obtained with the pump laser parameters and rotational temperature 
suggested by Boutu et al [ist. For the perpendicularly aligned molecules, the relative phase 
from both ensembles are almost identical and nearly independent of harmonic orders. This is 
in good agreement with the experiment. For the parallel case, our result with the pump laser 
parameters suggested in shows the phase jump starting near H31, which only mimics 
the experimental data. However, the phase jump slightly shifted to near H27 with the better 
aligned ensemble (solid line), bringing the result closer to the experiment. This indicates 
that the degree of alignment can play a critical role in determination of the precise position 
of the phase jump. A possible reason for the discrepancy is that the experimental setup 
was chosen such that the short trajectory is well phase-matched, whereas our simulations 
are carried out at the single molecule level with contributions from both long and short 
trajectories. In order to understand the origin of the phase jump and its dependence on 
degree of alignment, it is constructive to analyze transition dipole phase as a function of 
harmonic order for fixed alignments. This is shown in Fig. 12(b) for angular range from 
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20° to 60°. Clearly, the dipole phase shows a phase jump which shifts to higher order 
with larger angle. With a better alignment more contribution comes from small angles 
and the phase jump shifts toward smaller harmonic orders. Our calculations are carried 
out with a relatively low probe laser intensity of 1.25 x 10^^ W/cm^, the same as used in 



Boutu et a/ experiments [18|], in order to keep ground state depletion insignificant. As probe 
intensity increases our results indicate that the phase jump slightly shifts toward higher 
harmonic orders. Note that the phase jump cannot be reproduced within the SFA, where 
the insignificant phase difference between CO2 and Kr is caused only by the small difference 
in their ionization potentials (13.77 eV vs 14 eV), independently of alignment. 

V. VALIDITY OF THE TWO-CENTER INTERFERENCE MODEL 

Within the QRS theory, the structure of the HHG spectra directly relates to photo- 
recombination cross section. In particular, minima in photo-recombination cross section 
immediately result in minima in HHG spectra. It is therefore interesting to compare the 
positions of the minima in molecular systems under consideration with the prediction of the 



12|. According to this simple model, the minima 



simple two-emitter model by Lein et al 
satisfy the relation 

Rcos6 = {n+l/2)X'^^, n = 0,1,2,..., symmetric wavefunction, (66) 
= n\^^^^ ?T, = 1,2, anti-symmetric wavefunction (67) 

where X^^f is the "effective" wavelength of the continuum electron defined such that the 
"effective" wave vector is k'^^^ = -\/2?2 , with VL being the energy of the emitted photon. In 



other words, the energy is shifted by Ip with respect to the usual relation k = \j2{Vt — Ip). 

In Fig. 13 we show the projected internuclear distance RcosO vs electron "effective" 
wavelength at the minima of the photoionization differential cross sections for CO2, O2, N2, 
and H2. For completeness, we also plot here the results for H2' with different internuclear 



separations, reported in Ref . 37| . Remarkably, the minima from CO2 follow the two-emitter 



model very well, shown as the dotted line for the case of anti-symmetric wavefunction (due 



to the TTg symmetry of the HOMO). This fact 



shown in Ref. Il9| . For O2, which also has the rcg symmetry for the HOMO, the 



HHG measurement from aligned CO2 Q, [l5|, and more conclusive evidences have been 



las been first observed experimentally in the 
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FIG. 13: (Color online) Projected internuclear distance vs electron "effective" wavelength at the 
minima of the photoionization differential cross sections for CO2, O2, N2, and H2. Results for 
Hg" with different internuclear separations, reported in Ref. [s^], are also plotted for comparison. 
Predictions of the two emitter model are shown as dashed and dotted lines for symmetric and 
anti-symmetric wavefunctions, respectively. 

minima start to show up only at quite high energies (or shorter wavelengths). This is not 
surprising since the internuclear distance between the two oxygens is about two times shorter, 
as compared to CO2. The case of H2 also agrees reasonably well with the two-emitter model. 
However, N2 minima do not follow any simple pattern. This fact has been noticed earlier 



by Zimmermann et al 



97(1 . This clearly indicates that the two-center interference model is 



not guarantied to work a priori. 

Finally we show that the minima in the photoionization cross section of O2 can be ob- 
served in HHG experiments. Since the minima occur at quite high energies, it is better to 
use a driving laser with long wavelength of 1600 nm [93], as the cutoff will be extended to 
much higher energies without the need of using high laser intensities. In Fig. 14 we plot the 
photoionization cross sections for harmonic orders of 129, 159, 189 and 209, as functions 
of alignment angle. Note that the lower orders are already given in Fig. 5 (but in units of 
harmonics for 800 nm). The cross sections show a clear minimum for H159 near 37°, which 
slowly moves to larger angles for higher harmonics. This picture resembles the behavior for 
the CO2 case, but at much higher energies. Similar to CO2, a pump-probe scheme can be 
used to investigate the time delay behavior of the HHG yield. 

In Fig. 15 we show a typical behavior of the HHG yield as function of delay time between 
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FIG. 14: (Color online) Diffrential photoionization cross sections for O2 as functions of the align- 
ment angle. The photon energies are expressed in units of harmonic orders for 1600-nm laser. 




FIG. 15: (Color online) HHG yield from O2 for different harmonics as function of pump-probe delay 
time near quarter-revival. (sin^(20)) and (cos^(0)) are also plotted for reference (a). Calculations 
are performed with a 1600-nm laser pulse with peak intensity of 2 x lO^'' W/cm^ and 20 fs duration. 

pump and probe laser pulses for H49 and H69 in the lower plateau, H129 and H159 in the 
middle plateau, and H189 and H209 near the cutoff. The calculations were carried out with 
a 1600-nm laser pulse with peak intensity of 2 x 10^^ W/cm^ and 20 fs duration (FWHM). 
For the pump we use 800-nm pulse with 4 x 10^^ W/cm^ and 120 fs duration. We only focus 
on the results near the quarter-revival (near 3 ps), when the molecule can be most strongly 
aligned. This is in contrast to CO2 where most strong alignment is achieved near 3/4- 
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revival. The difference is caused by the different symmetries of the total electronic ground 
state wavef unctions, as discussed in Sec. 11(F). Clearly the behaviors are quite different for 
these energy ranges. In the lower plateau, the HHG yield follows (sin^(26')). This behavior 
is the same as for the case of 800-nm laser. In the middle plateau, the behavior starts to 
look similar to that of CO2, which follows the inverted of (cos^(^)). This reflects the fact 
that the photoionization DCS peaks have shifted to larger angles for higher photon energy. 
Finally, for even higher harmonics an additional peak occurs right at the quarter-revival. 
This is also similar to the case of CO2 near H31 for 800-nm, where the transition dipole 
goes through a minimum. Note that a quite high degree of alignment is needed in order to 
observe these additional peaks, which is not as pronounced as in CO2 near 3/4-revival. 



VI. SUMMARY AND OUTLOOK 



In this paper we give detailed description of the quantitative rescattering theory (QRS) 
applied for high-order harmonic generation (HH G) by intense laser pulses. Although HHG 
has been described by the three-step model 20|] in its many versions, including the Lewen- 



stein model 



22| . quantum orbits theory [99] since the mid 1990's, the calculated HHG spectra 



are known to be inaccurate, especially in the lower plateau region, as compared to accurate 
results from the TDSE. At the same time, numerical solution of the TDSE for molecules in 
intense laser pulses still remains a formidable challenge and has been carried out only for the 
simplest molecular system H2' so far. The QRS has been shown to provide a simple method 
for calculating accurate HHG spectra generated by atoms and molecules. The essence of 
the QRS is contained in Eq. fl54|) which states that harmonic dipole can be presented as 
a product of the returning electron wave packet and exact photo-recombination dipole for 
transition from a laser-free continuum state back to the initial bound state. The validity of 
the QRS has been carefully examined by checking against accurate results for both harmonic 
magnitude and phase from the solution of the TDSE for atomic targets within the single 
active electron approximation. For molecular targets, the QRS has been mostly tested for 
self-consistency in this paper. More careful tests have been reported earlier in Ref. for 
H^, for which accurate TDSE results are available. The results from the QRS for molecules 
are in very good agreement with available experimental data, whereas the SFA results are 
only qualitative at best. 
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Here are a number of the most important results of the QRS. 

(i) The wave packet is largely independent of the target and therefore can be obtained 
from a reference atom, for which numerical solution of the TDSE can be carried out, if 
needed. It can also be calculated from the SEA. The shape of the wave packet depends 
on the laser parameters only, but its magnitude also depends on the target through the 
ionization probability for electron emission along laser polarization direction. 

(ii) By using the factorization of HHG and the independence of the wave packet of the 
target structure, accurate photo-recombination transition dipole and phase can be retrieved 
:'rom HHG experiments. This has been demonstrated for atomic targets, both theoretically 



36|, |38| and experimentally [86||, and more recently for CO2 |8l| ]. 

(iii) Since the QRS is almost as simple as the SEA, it can be useful for realistic simulations 
of experiments, where macroscopic propagation needs to be carried out. Eor such simulations 
contributions from hundreds of intensities are needed. Existing macroscopic progagation 
simulations have rarely been done beyond the SEA model even for atomic targets. In the 
QRS the most time-consuming part is the calculation of the photoionizaion transition dipole. 
However, this needs to be done only once for each system, independently of laser parameters. 
The use of the QRS would significantly improve the SEA macroscopic results, but with the 
same computational effort as for the SEA. 

(iv) Because of the inherent factorization, the QRS allows one to improve separately 
the quality of the wave packet and of the transition dipole. In particular, one can include 
many-electron effect in the photoionization process, which has been performed routinely in 
the atomic and molecular photoionization research. 

The are several limitations of the QRS. Clearly, the method is not expected to work well in 
the multiphoton regime since QRS is based on the rescattering physics. Wave packets from 
the SEA becomes less accurate for lower plateau. We have seen that this has affected both 
harmonic yield and phase. This is probably due to the fact that the electron-core interaction 
has been neglected in the SEA during the electron propagation in the continuum. Remedy 
br this effect has been suggested by including some correction to the semi-classical action 



231] . We emphasize that the QRS so far only improves the last step of the three-step model 



by using exact photo-recombination transition dipole. Lastly we mention that similarity 
of the wave packets from different targets only holds for low to moderate intensities. Near 
saturation intensities, when depletion effect is large, the wave packets from targets with 
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different ionization rates would be very different. Nevertheless, one can still use the wave 
packet from the SFA, provided the depletion is properly included. 

The evidence presented in this work strongly supports that the QRS provides as a powerful 
method to obtain accurate HHG yield and phase from molecules under intense infrared 
laser pulses. As the technology of HHG generation improves, one can count on the retrieval 
of accurate photoionization cross section and phase from such measurement from aligned 
molecules. By using a pump beam to initiate a chemical reaction and a probe laser beam to 
generate high-order harmonics in a standard pump-probe scheme, the QRS would allow the 
retrieval of the transition dipole magnitude and phase over a broad range of photon energies, 
thus paving the way for extracting the structure information of the transient molecule, to 
achieve ultrafast chemical imaging with femtosecond temporal resolutions. In the meanwhile, 
with accurate single atom or single molecule induced dipole moments readily calculated, it 
is also possible to examine the effect of macroscopic propagation of the harmonics in the 
gaseous medium. 
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